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The observed neutron star mass distribution as a probe of the 
supernova explosion mechanism 
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ABSTRACT 

The observed distribution of neutron star (NS) masses reflects the physics of core-collapse su- 
pernova explosions and the structure of the massive stars that produce them at the end of their 
evolution. We present a Bayesian analysis that directly compares the NS mass distribution 
observed in double NS systems to theoretical models of NS formation. We find that models 
with standard binary mass ratio distributions are strongly preferred over independently pick- 
ing the masses from the initial mass function, although the strength of the inference depends 
on whether current assumptions for identifying the remnants of the primary and secondary 
stars are correct. Second, NS formation models with no mass fallback are favored because 
they reduce the dispersion in NS masses. The double NS system masses thus directly point 
to the mass coordinate where the supernova explosion was initiated, making them an excel- 
lent probe of the supernova explosion mechanism. If we assume no fallback and simply vary 
the mass coordinate separating the remnant and the supernova ejecta, we find that for solar 
metallicity stars the explosion most likely develops at the edge of the iron core at a specific 
entropy of S/Na ~ 2.8 fee. The primary limitations of our study are the poor knowledge of 
the supernova explosion mechanism and the lack of broad range of SN model explosions of 
LMC to solar metallicity. 
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1 INTRODUCTION 

At the end of the life of a massive star, the degenerate core made 
of heavy elements grows by episodes of nuclear shell burning until 
the Chandrasekhar instability initiates collapse. When the central 
matter reaches nuclear densities, the equation of state stiffens dra- 
matically as a result of the hard core repulsion of the strong force, 
and a shock wave is launched back into the supersonically collaps- 
ing outer mantle. The shock starts at an approximately constant 
mass coordinate of about 0.5 to 0.7 Mq, which depends mainly 
on weak-interaction physics and nuclear interactions, and little on 
the progenitor structure (e.g. iGoldreich & Weberll98oUYahilll983l ; 
Beth el 199d) . However, the mass of the iron core is too big to allow 
for a prompt explosion (e.g. lBethell990l) and the initial shock wave 
stalls due to neutrino losses and photodisentegration of heavy nu- 
clei and turns into a standing accr etion shock. This i s a robust fea- 
ture of supernova simulations (see ljanka et al.ll2007l for a review). 
The mass of the hot central object, the proto-neutron star (PNS), 
continues to grow by accretion. The PNS cools by emission of neu- 
trinos that are partially absorbed below the shock and may revive 
its outward movement to produce a core-collapse supernova (e.g. 



IColgate & Whitdl 19661 : fBethe & WilsorJI 19851) . The dense central 
remnant evolves either into a neutron star (NS), or into a black hole 



(BH) , if the explosion fai ls (e.g. Burrows 1986; Liebendorfer et al 



2001 



2011) 



Heger et al.1 120031 : iKochanek et alj 120081 : lO'Connor & OtJ 



if mass fallback during the explosion increases the mass 
above the maximum allowed NS mass (e.g. IWooslev & Weaveil 
1 19951 : 1 Zhang et al. 2008) or i f a phase transition occurs in the cool- 
ing NS (e.g lBrown & Bethelll994l;lKeil & Jankall 19951) . 

As the supernova explosion develops within the inner core 
of the star, the overlying layers of stellar material prevent di- 
rect observations of the region aside from neut rinos or gravity 
waves emitted by very nearby superno vae (e.g. lOtt et al.l 120041 : 
lYtiksel & Beacomll2007r . lOtt et alj|2012l) . Thus, our understanding 
of the relevant processes is largely based on numerical simulations, 
which generally do not produce explosions for prog e nitors more 
massiv e than about 15 Mp (e .g. Kitauraetal. 2006; Buras et al 
2006b:; iMarek & Jankj 120091: ISuwa et al.l l201Ct iTakiwaki et al 



120121 ; iMuller et alj |2012|T Model explosions fail because the 



neutrino luminosity never reaches the critical value necessary 
to turn the quasi-stationary accretion shock into an outgoing 
explosion in a non-rotating progenitor (e.g. iBurrows "&~ Gosh v 
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1993l;lYamasaki & Yamadall200ll200d:lMurphY & Burrowsll2008l; 



Nordh aus et al.l20ld : |Peicha & Thompsonl2012l : lFernandezl2012h . 

These models with fail ed explosions continue to accrete until the 
forma tion of a BH (e.g. lLiebendorfer et alJkoOlUO'Connor & Ottl 
J). This "supernova explosion problem" has remained unsolved 
for more than four decades. Present day studies are focused on 
whether mu l ti-dim e nsiona l effects will aid the explosion (e.g. 
Herant et alj 1 1991 fl994l: iBurrows etaLl Il995t Jjanka & MulQ 

Warrenl 120021 12004 
Ijjf iMurphv & Burrows! 



19961: iFrver & Hegeri |200(A lFrver& 
Ohnis hi et al. 120061: llwakami et alj " 



200a iMarek & J anka 2009]; iNordhaus et all 120101: lHanke et al.l 
201 lLlTakiwaki et al.l2012l) . As these simulations address primarily 



non-rotating progenitors, it is worth mentioning that the combina- 
tion of sufficiently rapid rotation and strong magnetic fields can 
potentially explod e a wider range of prog e nitor stellar masses (e.g . 
Svmbalistvl 1 1984 l Akivama et alj 120031; IThompson et alj 120051 : 
Burro ws et al ■l2007l:lDessart et a l. 2008, 2012h. 



In addition to the mass at core bounce and the mass ac- 
creted through the shock before explosion, the final mass of the 
core-collapse remnant is determined by the amount of fallback 
during the explosiorQ. Because self-consistent supernova simula- 
tions do not genetically explode, the usual approach for study- 
ing fallback, supernova nucleosynthesis, and their observational 
consequences is to initiate an artificial explosion by depositing 
enough momentum or energy at a specific mass coordinate of the 
progenitor to match the final supernova energy or nickel yield 
(Wooslev & Weaver 1995; Thielemann et al. 1996; Timmes et al. 
1 19961 : IZhang et alj|2008l : iDessart et ai]|201ll) . The boundary con- 
ditions are adjusted for each progenitor star to give an explosion 
with the desired properties. The amount of fallb ack in these mod- 
els depends on the boundary conditions ( MacF adven et alj|200ll) 
and the structure of the progenitor. Broadly speaking, these simpli- 
fied models predict that low-metallicity stars explode as supernovae 
when they are still hot and blue. Their compact envelopes thus gen- 
erate s tronger reverse shocks, more fallb ack, and higher mass rem- 
nants dChevaliej [T989; Zhan g et alj|2008l) . Due to mass loss, solar- 
metallicity stars have less massive final envelopes, which generi- 
cally lead to less fallback, and the stron gest reverse shocks are pro - 
duced at the helium/hydrogen interface dWooslev & We aver 1995). 
Stars without hydrogen envelopes, such as Wolf-Rayet stars, expe- 
rience little fallback and thus always produce NSs if the supernova 
mechanism is successful. In any model, more energeti c explosions 
produ ce less fallback and smaller remnant masses (Zhang et al. 

2008) . In addition to increasing the mass of the remnant, variations 
in fallback could intro duce an element of stochasticity into the NS 
or BH mass function dOzel et al.ll2012l) . 

Observational constraints on the causal chain that links mas- 
sive stars, the supernova mechanism and their remnants are diffi- 
cult to obtain and have many open problems. Examinations of pre- 
explosion images of type Hp supernova^ suggest tha t the progen- 
itors have initial masses lower than 16.5 ± 1.5 Mq dSmartt et al.1 

2009) and , more generally, there is a dearth of high-mass progen- 
itor stars dKochanek et alj|2008l) . The upper limit on the Hp pro- 
genitor mass is surprising, becau se red supergiants wi th masses 
of up to 25 Mq are observed dLevesque et al.1 120051) and are 
thought to explode and produce NSs dHeger et al. 1 120031) . While 



1 The fraction of the PNS m ass lost due to neu trino-driven wind in the 
explosion is negligible (e.g. IThompson et al . 2001). 

2 The plateau in the light curve indicates a presence of thick hydroge n 
envelope and a red supergiant progenitor (e.g. Chevalier 1976; Amett 1980). 



the possibility that the se stars do not explode at all is intriguing 
(Koch anek et alj [20081 ). the stati stical significance o f this "red su- 
pergiant problem" is only 2.4a dSmartt et al .1120090 . Other expla- 
nations involve red supergiants exploding as other types of su- 
pernovae due to differences in mass-loss ra t es dSmith et alj|2009l : 
lYoon & Cantielldl20ld: [Moriya et al.ll201 ll; lGeorgy||2012l) . or bi- 
nary evolution (Eldridge et a l.ll201 lHSmith et al.ll201lF l 



As the remnant mass function encodes information about 
the structure of the progenitors and the supernova explosion 
mechanism, considerab le attention has been devoted to measur- 
ing NS and BH mas s es {Finnlll994[Thorsett & Chakrab; 



20 id ; 



Valentim et 



f !999 
2011 



Schwa b et al] l20ld: IK iziltan et 
Zhang et alj|201 ll ; IQzel et alj|2012h . The most precise mass mea 
surements come from binary systems where one component is a 
pulsar and there are accurate measurements of at least two post- 
Keplerian parameters. A special group of such systems are double 
NS binaries (DNS), where the first-born NS is observed as a recy- 
cled pulsar with a rotation period of 20 ms < P < 200 ms, and the 
companion NS is the result of a supernova from the initially less 
massive secondary star (e.g . Bhattacharva & van den Heuvel 1991; 
Porteg ies Zwart & Yun gelson 1998; B urgav et al.ll2003 ; Lvne et all 



2004) . Six such systems are currently known, giving 12 pre- 
cise mass measurements. The masses of these NS s cluster at ~ 
1.35 Mp) with a disper sion of only ~ 0.06 Mq dKiziltan et alj 
l20ld : IQzel et alj |2012|) . and the mean masses of the pulsars 
and their companions differ by only 0.03 Mq d6zel et aljHoTl) . 
Present theories argue for very little mass transfer in these systems 
once the first NS is formed, so both NS masses basically represent 
the birth masses of the NSs. The mean observed masses are sig- 
nificantly higher than the Chandrasekhar mass of the pre-collapse 
core, which may i ndicate growth by fallback during the supernova 
dKiziltan et alj|201fj) . The tight mass distribution, however, argues 
against significant fallback, and the properties of the distribution 
have also been attributed to the particular evolution history tha t 
leads to their formation d6zel et al.1 12012|) . ISchwab et ail j201Ch 
proposed that a third of the DNSs formed as a result of electron- 
capture supernovae as evidenced by their lower masses, with the 
remaining higher mass systems forming as a result of a Fe-core 
collapse. 



IKiziltan et all feOld) and lQzel et ail d2012l) have carefully fit 
simple analytic models to the DNS mass distributions. Here, we 
take the additional step of using the mass distribution of NSs to 
probe the physics of supernova explosions. We replace the para- 
met ric models of obse rved NS masses used bv lKiziltan et al. |{20l3) 
and lQzel et alj d2012h with predictions of NS masses based on the 
actual physics of progenitors and supernova explosions so that we 
can directly constrain, compare and assess the validity of physi- 
cal models for the explosion. We present a Bayesian formalism 
that quantitatively compares different predictions for the remnant 
mass function to the DNS data. We compare different DNS pro- 
duction mod els using the artificial supernova explosion models of 
IZhang et al.1 d2008l) with different explosion energies and differ- 
ent progenitor metallicities. In Section [2] we outline our Bayesian 
framework for comparing the NS production models with the ob- 
served data. In Section [3] we present our results and outline ex- 
tensions of our model that may further constrain the underlying 
physics. In Section[4] we discuss our results and their implications 
for the supernova explosion mechanism. 
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2 STATISTICAL MODEL 

In this Section we outline a general Bayesian statistical model to 
quantitatively evaluate different hypotheses about the origin of the 
DNS mass distribution. We choose a Bayesian framework because 
it easily allows for a simultaneous comparison of multiple mod- 
els with different numbers of parameters, yields best-fit parameter 
estimates, and naturally incoiporates prior knowledge. We start by 
outlining a procedure for calculating the posterior probability dis- 
tributions and then we formulate several hypotheses to be evalu- 
ated. We also describe the data and the underlying physical model. 
A major limitation for these models is the very limited availability 
of supernova explosion models - even those employing a simple 
piston at a fixed composition jump or mass cut - as a function of 
mass and metallicity. 



2.1 General considerations 

According to Bayes theorem, the posterior probability of hypothe- 
sis H with internal parameters given data D, P(0H\D), is equal 
to the prior probability P(8H) multiplied by the marginal likeli- 
hood P(D\8H ) that D arose from hypothesis H, 



P(8H\D) oc P(8H)P(D\8H). 



(1) 



If the data D are composed of N individual measurements and 
the i-th measurement is characterized by a probability density in 
an observed pair of masses M = (Mi, M2), Pj(D|M), then the 
marginal likelihood of hypothesis H is 



P(D|0P) = Y[ / Pi(D\M)Pi(M\0H)dM, 



(2) 



where Pi(M|0P) is the probability that the given value of M 
occurs for the parameters of H for the i-th measurement. For 
a given hypothesis H, different values of the parameters yield 
different probabilities of the data P(D\OH) such that we can de- 
termine the "best-fit" parameters and their confidence intervals 
based on the posterior probability distribution P(8H\T>). 

Suppose that we have two hypotheses Hi and H 2 parameter- 
ized by their individual parameter sets 0i and 8%. Which of the 
two hypotheses better describes the data? Within the framework of 
Bayesian analysis, the relative "probability" of the two hypotheses 
is given by the Bayes factor 



-B12 = — = 



B 1 _ /P(D|0iifi)P(0i#i)d0i 



B 2 J P{D\8 2 H 2 )P(8 2 H 2 )d8 2 ' 



(3) 



Note that only the ratio of E>i and B2, B12, has any meaning and 
that it can be extended to an arbitrary number of hypotheses. Proper 
calculation of B12 also requires that the individual probabilities in 
Equation l(2j are properly normalized with J Pj(D|M) dM = 1, 
/P(M|0P")dM = 1, and fP( 8H) d0 = 1 over the relevant 
ranges of M and 0. |jeffrevsl l l 19831) groups values of P12 in several 
categories: B12 > W 1 ^ 2 implies that hypothesis Hi is "substan- 
tially" better than H 2 . If B12 > 1 2 , then the ev idence against 
H2 and in favor of Hi is decisive. Jeffreys] dl983l) also gives ta- 
bles to approximately relate P12 as a function of number of the 
parameters in to a more commonly used \ 2 difference, specifi- 
cally B12 oc exp(-Ax 2 /2). 

The Bayesia n statistical mode l we pr esent here is similar to the 
one developed bv lOzel et al.l fe0ld . l2012h with a key difference: in- 
stead of using a phenomenological description based on a paramet- 
ric function (in their case, a Gaussian), we will tie the observed NS 



masses directly to physical calculations of remnant masses based 
on supernova physics and the progenitor structure. This allows us 
to quantitatively compare different scenarios for the origin of the 
NS mass distribution. 



2.2 NSs as members of a binary 

The masses of the two NSs in a binary system are not independent 
and reflect the binary initial mass distribution, any mass transfer 
processes that occurred during the system evolution, and the super- 
nova physics. We assume that one observation yields a pair of NS 
masses of the binary, M = (Mi, M2), where Mi is the mass of 
the recycled pulsar, and M2 is the mass of the companion. In Equa- 
tion {2}, Pj(D|M) is the probability of observing the i-th pair of 



P(D|M) = JV(Mi,Mi 1 i,a 1> i)Af(M 2 , M 2 ,i,a 2 ,i 



(4) 



and (Mi,i,M2,»), &2,i) are the measured masses and their 

uncertainties for the i-th DNS system. Here, N are Gaussians de- 
fined as 



Af(x,n, a) 



exp 



2a 2 



(5) 



V2na 2 

Equation I0 assumes that there are no correlations between the 
two NS mass measurements in the binary, although these could be 
included. More complicated models of the probability densities of 
the observed masses can be included as well. 

An important issue for the calculation of P;(M|0P") is the 
assignment of the DNS components to the original primary and 
secondary stars in the binary. There are two mutually incompatible 
possibilities, specifically, either the recycled pulsar came from the 
primary sta r and the co mpanion from the secondary, or the reverse. 
Following IPressI (1997), we account for these two probabilities by 
expressing Pi(M|^# max ) as a combination of these two mutually 
incompatible hypotheses 



(6) 



Pi(M\8H) = J P( Pl )[ Pl 3 g (M 1 ,M 2 \8H) + 
+(l- P i)&>(M 2> Mi\8H)] d Pl , 



where pi is the probability that the recycled pulsar in system 
i came from the primary star and P(pi) is the prior on pi. 
£?{Ma , Mb I OH) is the probability density distribution of the pair 
of remnant masses (Ma, Mb) where Ma corresponds to the rem- 
nant mass of the primary star and Mb to the that of the sec- 
ondary star. For an uniform prior P{pi), the final probability is 
a simple average of the two options. However, the binary evolu- 
tion models generally require that the recycled pulsar originated 
from the primary star and hence the prior P(pi) is strongly peaked 
at pi = 1. Based on the binary evolution models, we choose 
P(pi) = S(pi — 1) for all systems. In Section [3~T1 we investigate 
the appropriate choice of P(p%) for the individual systems. 

We consider two forms of the probability J 2 (Ma, Mb\OH). 
First, as a counterpoint to the more complicated models, we con- 
sider an "independent" star model where each star is independently 
drawn from the IMF. This simple model is conceptually similar to 
the parametric mod els that assume no co rrela tion between the s tars 
in a binary used by Kizilt an et ail 120 l(t) and lQzel etal] d2012l) . In 
this model the primary and secondary probability distributions are 
independent, 9^(M A ,M B \8H) = 3^(Ma\0H)&>(M b \0H), 
and ^(AIIOH) is the probability distribution of remnant masses 
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of a single star 



w max 

&>(M\OH)ix J P{M)M[M, M'{jV),a thKa ] AJK. (7) 



Here, we assume that NSs are produced by stars with initial 
massefl JK between ./idm and ^# ma x, where JK m i n is fixed 
and JIijx is a parameter (^max £ 9)- P(^t) is the probabil- 
ity of progenitor mass ./£ , which we assu me to be a powe r law, 
P{J£) oc J(~ a , with a = 2.35 to match lSalpeterl jl955t) . The 
function M {./£) provides the remnant mass M for the given pro- 
genitor mass ^£ (see Section l2?3t . The independent model is sym- 
metric and thus Equation l[6j will give the same marginal likelihood 
with no dependence on P(pi)- 

Second, we consider a genuine binary distribution that is de- 
fined as 

^max -di A 

&>(M A , Mb \0H) oc J cUCt J AJ^ b P{J(a)^^- X 

xM[M A , M'{JK A ), <T theo ]Af[M B ,M'(^ B ), a thco ], (8) 

where we assume that the primaries are drawn from a Salpeter 
IMF, P(./#a) oc ,/^" 2 ' 35 , and the secondaries are drawn from a 
distribution P(q) of mass ratios q = .MbI-^a- We assume ei- 
ther uniform P(q) for 0.02 $; q ^ 1.0, or a population of "twin" 
binaries with half of binaries distributed uniformly in the interval 
0.9 ^ q ^ 1 and the other half uniformly distributed for 0.02 ^ 
a < 0.9 jPinsonneault & StanektoOfj ; iKobulnickv & Frvej|2007l ; 
lKochanelJ T2009). We normalize the mass ratio distribution as 
/o 1 2 P{<i)dq = 1 and we drop systems with secondaries with 
< ^min that would produce NS-WD binaries. Systems with 
more massive primaries thus produce a higher relative fraction of 
DNSs. We neglect all binary evolution processes that could mod- 
ify the relation between the initial and remnant masses M ' 
because the mass transfer in a DNS progenitor binary system oc- 
curs after the main sequence evol ution, which fixes the size of the 
helium core of the primary (e.g. Bhattacharva & van den Heuvell 



Il99ll : IPortegies Zwart & Yungelsonlll998l) . We thus assume that 
M'(y$) is the same for primaries and secondaries. Again, ./# max 
is a free parameter (^ max € 0). 

We chose the form of Equations <[7j — <[8j f° r several reasons. 
First, the function M'{j$) is usually tabulated only for a dis- 
crete set of Jt( and we need &>(M\0H) and 0»{Ma,Mb\9H) 
to be continuous and smooth. This is because the likelihoods of 
many of the NS mass measurements are sharply peaked and we 
do not want our results to be sensitive to the exact position of 
the NS mass with respect to the discrete mass models of the the- 
oretical studies. Second, the width of the kernel, athco, can be 
interpreted as the uncertainty in the theoretical NS masses either 
due to progenitor structure, NS growth during the accretion phase, 
or stochasticity in the amount of fallback. In principle, one can 
have (Tthco = c"thco(.^) and make the NS mass uncertainty de- 
pend on the progenitor mass. If crtheo is too small, £P(AI\9H) 
and £P(Ma, Mb\9H) will have many individual peaks, while if 
it is too large, the structure in the NS distribution will be smeared 
out. We varied 0.01 Mq ^ (Tthco ^ 0.05 Mq and found that for 
higher values of er t h C o the relative probabilities of the models were 



3 Throughout this paper we denote progenitor masses as jM and remnant 
masses as M. 



smaller. However, the ordering of the models did not change. We 
choose uthco = 0.025 Mq as a rather arbitrary compromise be- 
tween the two extremes. This width is several time s smaller than 
the typ ical width of the DN S mass distributions of Kizilt an et al.l 
feOlCh and lOzeletai]j2012h . 

The last quantity we need to evaluate Equation QJ is the prior 
on Xax, f (^Cnax). Since Jmax attains only positive values and 
we do not have any physical constraints, we set the prior to be uni- 
form in ln./# max for 10 Mq ^ Xa, ^ 100 Mq, where the up- 
per limit corresponds to the approximate maximum mass of a star. 

2.3 Data, underlying models and implementation 

In order to calculate P(6H\~D) for the independent and bi- 
nary models, we need the mapping between the initial progenitor 
mass and the fina l remnant mass M'(^M\ We use the results of 
IZhang et al.l 120081) summarized in Table [T] who obtained NS and 
BH mass distributions for primordial (Z = 0) and solar metal - 
licity {Z = Zq) progenitors by positioning a piston at a partic- 
ular mass coordinate and injecting enough momentum to obtain 
an explosion with the desired ejecta kinetic energy E at infinity. 
The pistons were positioned either at the point where the entropy 
S/Na = 4 &b, which corresponds approximately to the base of the 
oxygen burning shell, or at the edge of the deleptonized core ("Y e 
core"), which is located deeper in the star where the electron frac- 
tion Y e decreases due to electron captures on protons. This radius 
roughly corresponds to the iron core. We also consider remnant 
masses that correspond to the Y c core and S/Na = 4 fee masses 
with no fallback. Looking at Table [T] we see that the model cal- 
culations do not extend all the way to the minimum mass for su- 
pernova explosion „# m in. It is expected that fallback is negligible 
for these low mass sta rs and that the remn ant mass is equal to the 
core mass. Following IZhang et al.l J2008h . we extend the proper- 
ties of the 10 Mq stars down to ^# m in for primordial composition 
stars. For solar metallicity and the piston at S/Na = 4 fee, we set 
M = 1.37 Mq for 11 < J£ < 12 M and M = 1.35 Mq for 
9.1 ^ . // < 11 Mq. For solar metallicity and the piston at the Y c 
core, we set M = 1.32 Mq for 9.1 < Jt < 12 Mq. All remnant 
masses were corrected for the loss of binding energy i?bmd due to 
neutrino emission during the supernova event using the approxima- 
tion 



Pbin 



0.075 M, 







M s 



Mr 



(9) 



where M grsv is the gravitational mass o f the remnant after the cor- 
rection for -EWnd dTimrnes et al.lll99r3) . The co mmonly assumed 
value of 0.084 Mq for the leading factor from lLattimer & Yahill 
( 1989) diff ers slightly from Equation (|9). We adopt the sample of 
DNS from lOzeletafl J2012h . which is reproduced in Table [2] for 
convenience. The sample consists of 6 NS binaries, that yield 12 
precise NS mass measurements. 

Because we used Gaussians for P(DjM) and the kernels ap- 
pearing in Equations (0 and l[8}, we can evaluate the integrals over 
M in Equation ([2} analytically by swapping the order of integra- 
tion. This greatly speeds up the calculation, especially for the bi- 
nary models. Integrals over ^£ in Equations 10 and ([8} were eval- 
uated using the midpoint rule centered on the given progenitor j$ 
and with equal to half the distance in ./£ to the nearest progen- 
itor models. We did not use a more sophisticated integration method 
because the NS mass distribution is not intrinsically smooth and 
because there are significant jumps in M' between progenitors of 
similar mass. We use 5 to 10 points for the mass ranges where no di- 
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rectly calculated progenitors are available (9.5 ^ ^ ^ 10 M© for 
Z = and 9.1 sC Jt < 12 Mq for Z = Z q), We assume that th e 
maximum gravitational NS mass is 2.0 M© dDemorest et al.l2010h . 
and we thus do not include in the calculation of ^(M\0H) or 
^>(M A ,M B \0H) any progenitor producing M'{Jf) > 2.0 M . 
These progenitors are assumed to yield BHs. 



3 RESULTS 

We first discuss the DNS mass distribution and the ambiguities in 
associating a NS in a binary with a progenitor (Section [3T). Then 
we examine a range of hypotheses on the origin of the NS mass 
distribution (Section [3.2t . We also present several extensions and 
limitations to our analysis (Section [3.3t . 

3.1 Properties of the binary model 

In Figures [TJ and [2] we show examples of the probability distribu- 
tion of NS masses originating from a binary &{Ma, M B \^max) 
with a uniform P(q) along with the distributions marginalized over 
Ma or Mb - The model in Figure [TJhas the highest relative proba- 
bility of all models considered in Section POl For comparison, Fig- 
ure|2]shows a model with fallback that has relative probability B12 
(Eq. (3)) lower by a factor of ~ 10 2 S . Examining Equation {8]l, 
we see that the primary and secondary mass ranges producing NS 
are quite different. For uniform P(q), the probability distribution of 
primary masses jMa is proportional to ^^ a (l— ^Kmin/^A), and 
primaries with masses close to ./# m i n do not contribute to the dis- 
tribution of NS masses Ma because these systems mostly produce 
NS-WD binaries. The primary distribution producing NSs peaks at 
(q + l)./iCnin/ ot ~ 1.43^# m i n for a Salpeter IMF. The distribution 
of secondaries producing NSs is proportional to — 
for a flat P(q). Here, the cutoff is at high masses, while the dis- 
tribution of the lowest mass progenitors is almost Salpeter. These 
analytic estimates immediately show that NS masses originating 
from the primary and secondary stars of a binary represent differ- 
ent progenitor mass ranges. For example, for ^# m in = 9.1 Mq 
and ,-# milI = 25 Mq , the mean progenitor masses of the primary 
and secondary components are 16.5 and 12.8 M©. This explains 
why the marginal distribution of Ma in Figure Q] is not peaked at 
low Ma and has a stronger secondary peak at Ma ~ 1.6 M©, 
when compared to the distribution of Mb ■ The secondary peak is 
much higher for models that include fallback (Fig. |2). There are 
no observed DNSs with masses in this second peak, which leads to 
a preference for models with no fallback (Section [3.2l . The secon- 
daries ,/£ B have progenitor distributions close to the IMF. Thus, the 
highest peak for secondaries is at 1.22 Mq, which is the assumed 
gravitational NS mass for stars with 9.1 ^ ^ 12 Mq. 

We see from Figures[T]and[2]that most of the probability is in 
the region where the primary produce s a more massiv e NS. How- 
ever, since M'(-#) is not monotonic (Zhang et al. 2008), there is a 
small probability that the more massive NS originated in fact from 
the less massive progenitor. Generally, the farther the NS mass pair 
is from the diagonal (Ma = M B ), the smaller the probability that 
the more massive NS came from the less massive progenitor. In 
our formalism we can estimate whether the mass difference be- 
tween the two NSs is enough to distinguish between an NS orig- 
inating from the primary or the secondary, or essentially whether 
the prior P(pi) = — 1) in Equation $6^ is appropriate. Fig- 
ure [3] shows the ratio of probabilities that the millisecond pulsar 
originated from the more massive progenitor, 3^{M\, A^l^diax), 



as compared to the reverse ^(M2, Mi |^ max ). In Figure [3] we 
have marginalized over explosion energies and J mM . We see that 
if I Mi — M2I < 2<7 t hco, our model cannot distinguish between 
the primary/secondary origin of the millisecond pulsar based on 
the masses alone. The case of J1906+0746 is peculiar (because the 
pulsar is significantly less massive than the companion) and our 
results show that it is unlikely to have originated from the more 
massive progenitor. In agreement with lLorimer et al.l (120061) who 
give a very small characteristic pulsar age (see also our Table |2](, 
we propose that the observed pulsar in J 1906 +0746 comes from 
the less massive secondary star and we set the prior on pi in Equa- 
tion ® to be P(pi) = 8(pi) for this system. We show this alterna- 
tive assignment as an open circle and the dashed lines in Figure Q] 
This alternative assignment increases the relative probability of the 
binary models by a factor of ~ 3 to ~ 300. 

3.2 Comparison of the individual models 

Next we evaluate the relative probabilities of individual models. We 
specifically discuss the differences between the independent and bi- 
nary models, the explosion energy, and the position of the piston. 
We find that there is little difference in relative probability between 
the uniform and twin mass ratio distributions so we only discuss the 
uniform P(q) model, which has slightly higher probability for solar 
metallicity. By comparing the relative probabilities of models with 
free ^ ma x to models that include all progenitors (equivalent to set- 
tingP(.^ max ) = <5(>f max -lOOM ) inEq. Q), we also find that 
models with ./# max as a free parameter are not significantly pre- 
ferred. The inferred values of .-# max range from 14 M© to 35 M© 
depending on the method used to infer the "best-fit" value, but with 
confidence intervals covering most of the allowed range for ^ max . 
Here, we show models marginalized over ^# max , although models 
simply fixing ,/# max = 100 M© give essentially the same results. 
Finally, for the purposes of this Section, we do the calculations with 
the puls ar in J 1906+0746 a ttributed to the secondary star, consis- 
tent with lLorimer et aL l d2006l) and our discussion in Section [3~T1 

Figure|4]shows the relative probabilities (Bayes factors) of our 
models as a function of explosion energy. Focusing first on the pri- 
mordial composition models (left panel of Fig. [4), we see a general 
trend of increasing model probability for explosion energies up to 
about E ~ 2 x 10 51 ergs, after which the relative probability is 
essentially constant. The explanation is that higher E explosions 
produce less fallback and hence reduce the number of high-mass 
NSs that expand the overall NS mass range. This is confirmed by 
the models that include only the mass of the core with no fallback 
(horizontal lines in Figure |4) that match the probabilities of the 
high-i5 models. For solar compo sition (right pa nel of Fig. |4), the 
effect of E is not clear. Since lZhang et al.l (12008) investigated only 
two values of E, the probability ratios are not large (factor of ~ 3). 

In all cases, the relative probabilities for binary versus in- 
dependent mass models is between 5 and 200, which suggests a 
strong preference for binary models unless we assume that the pul- 
sar in J1906+0746 came from the initially more massive star. In 
this case, the relative probability of the binary models decreases to 
a factor of ~ 3 for Z = and for Z = Zq some of the mod- 
els even disfavor the binary models. This relative change was ex- 
pected based on Figure [3] If we assume uniform P(pi) in Equa- 
tion ®, the relative probability of binary models is again only a 
factor of ~ 3 higher than for the independent mass model. Binary 
models are significantly punished if there is a system with Ma sig- 
nificantly lighter than Mb if that is not allowed by the underlying 
remnant mass model. Correct treatment of the primary/secondary 
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millisecond pulsar mass M A [M Q ] ?(¥jll ma J 

Figure 1. Contours of probability £?(Ma , Mg |./#max) for having a DNS with masses Mj\ and Mg that originated from the more massive primary and less 
massive secondary stars of the original binary, respectively, for the most probable model discussed in Section |X21 The model has uniform distribution of q 
and the remaining parameters are given in the upper right corner of the plot. Contour labels indicate the total probability that lies outside the contour (i.e. the 
contour labeled as "0. 1%" encloses 99.9% of the total probability). Circles indicate pairs of DNS masses. For the red systems, the pulsar is always assumed 
to arise from the primary. The green circle is for J 1906 +0746 where we show both possible assignments (see text). The upper panel shows the distribution of 
remnant masses Ma from the primary, with the lines marking the observed values of Mi . For J1906+0746, we show both cases, where the pulsar came from 
the primary (solid) or the secondary (dashed). The right panel shows the distribution in Mg . 



assignment of the millisecond pulsar and companion is crucial for 
properly evaluating the relative probabilities of the independent and 
binary models. 

Figure |4]also shows clear differences in the relative probabil- 
ities of the different piston positions. For Z = 0, the pistons at 
the Y c core are strongly disfavored, because for low-mass progeni- 
tors the masses of the Y c cores are too low. For solar composition, 
the situation is reversed. Models with the piston at the Y c core are 
significantly more likely than those putting it at S/Na — 4 ks, be- 
cause they reduce the number of NSs with M > 1.5 Mq. Further- 
more, we see that the highest probability models essentially corre- 
spond to those with no fallback. There are small changes in the rel- 
ative probabilities between the two binary mass ratio distributions 
(P(q)) and the primary/secondary assignment for J1906+0746. 
But the best overall model is the one with no fallback, uniform 
P{q), and remnant masses equal to the Y c core mass. The two pis- 



ton positions used bv lZhang et all (2008) are somewhat arbitrary, so 
in Section|4]we address the question of whether some other piston 
position would produce better agreement with the observations. 

While in the rest of this paper we examine the relative proba- 
bility of different hypotheses for the DNS mass distribution in the 
Bayesian sense, at this point we make a "frequentist" diversion and 
compare the observations to the models in an absolute sense. We 
have argued that the individual NS masses in a DNS system reflect 
the binarity of the original system in the sense that the primary and 
secondary progenitor and NS mass distributions are different. Even 
though one of the components is observed as a pulsar, the ambigu- 
ity of which component originated from the primary or secondary 
star of the binary remains present for some systems unless more in- 
formation is added (J1906+0746). In order to compare the models 
globally, we construct in Figure [5] cumulative distributions of the 
DNS total mass (Mi + M 2 ) and the mass difference (|Mi - M 2 \), 
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Figure 2. Same as in Figure[T] but for a model with significant mass fallback. The model has uniform P(q) and the remaining parameters are given in the 
upper right corner of the plot. The relative probability B12 of this model is a factor of ~ 10 2 ' 8 lower than for the model in Figure[T](see Section [Ol for more 
details). 



along with the predictions from Zhan g et al .(2008) models coupled 
to the scenarios described in Section[2]with ^ max = 100 Mq. 

The distribution of mass differences in the left column of Fig- 
ure [5] shows that getting the observed mass differences is entirely 
within the range of the models, although most of them predict a 
much broader distribution of mass differences. The twin models 
(blue and green lines) generally give smaller mass differences than 
their uniform P(q) counteiparts. Note that the no fallback, Y c core, 
solar metallicity model with the highest Bayesian relative probabil- 
ity (thick yellow dashed line in the lower left panel) shows almost 
perfect agreement with the observations. We see in the cumulative 
distribution of total masses in the right column of Figure |5]that the 
observed distribution is much narrower than all theoretical predic- 
tions. The narrowest cumulative distributions are again produced 
by solar metallicity models with no fallback and remnant masses 
given by the Y e core masses (thick dashed lines). Implementing 
a cutoff for progenitor mass ./# max makes the distributions nar- 
rower. However, this is only weakly favored by the observations. 
Additionally, the observed minimum total mass of ~ 2.6 Mq is 
markedly higher than any of the minimum total masses predicted 



by the models (2.2 to 2.4 M©). The shift is 0.2 to 0.3 Mq, which 
is much higher than the mass necessary to recycle the pulsar to the 
observed spin periods. 

At this point in a "frequentist" analysis, we would compare the 
two cumulative distributions using a Kolmogorov-Smirnov test to 
ascertain whether the observations are compatible with the models 
in an absolute sense. However, mass differences and total masses 
are only particular aspects of the full 2D distributions. If we exam- 
ine FigureQ] which shows &>{M\ , M 2 1 l mal = 100 Mq ) for the 
model with the highest Bayesian relative probability, we see that 
all 6 binaries lie within the probability contour containing 75% of 
the probability (if we assume that the pulsar in J 1906+0746 came 
from the initially less massive progenitor). The chance of having no 
system outside this contour is 0.75 6 = 0.178 and we would typi- 
cally expect 4.5 ± 1.1 systems within this contour given 6 systems. 
This suggests that the highest Bayesian probability model repre- 
sents the data quite well - the fact that there is no DNS system 
with M > 1.5 Mq is likely only a statistical fluctuation. How- 
ever, none of the three additional DNS systems with accurate total 
masses (Fig. [5] right panels) has a total mass of about 2.4 Mq, 
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Figure 3. Probability that the pulsar in the DNS came from the initially more massive progenitor with respect to the reverse. DNS system names are given at 
the bottom while the individual NS masses are given at the top with the pulsar masses above the companion masses. The symbols indicate the median, while 
the error bars shows 1 and 2<r contours. The results are shown for the two piston positions, two metallicities and the two binary models (symbols are explained 
in the plot), and are marginalized over all values of E and ./# max , weighted by P(.»# max |D). Filled squares and circles mark the median, while the error bars 
show 1 and 2cr quantiles. 
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Figure 4. Relative probability of the different models for the origin of the DNS mass distribution as a function of the explosion energy and marginalized over 
^max- We show primordial composition (left panel) and solar metallicity (right panel) models for the piston positioned at S/N^ = 4 kg (solid lines with 
filled symbols) and at the Y e core (dashed lines with open symbols). For the sake of clarity, we show only the independent mass model and the uniform P(q) 
binary model. The horizontal lines indicate results for the models with NS masses equal to the mass at the piston position and no fallback. 



© 0000 RAS, MNRAS 000, 000-000 




Figure 5. Cumulative prob ability distributions of the mass differences (\Ma — Mg\, left panels) and the total masses (Ma + Mb, right panels) of the DNS 
systems. The results of the lzhang et alj fc008T) models with ^max — 100 Mq are shown with thin lines for Z — (top panels) and Z — Zq (middle 
panels). The bottom panels show solar metallicity models with no fallback. The observed cumulative distributions of \Mi — M2 1 and Mi + M2 are shown 
with the thic k black li ne. The distribution of total masses in th e right panel includes additionally DNS systems J 15 18+4904 ( Janssen et al. 2008), Jl 81 1 — 1736 
iCorongiu et al.l2007l) . and J1829+2456 (Champion et al. 200j that have accurate total (but not individual) masses. 



which suggests that the finer features of the DNS mass distribution 
might not be entirely reflected in the theoretical models. More DNS 
systems with accurate masses of both components are necessary to 
address this question. 

3.3 Extensions and limitations 

There are numerous possible extensions to the analyses presented 
in this paper. For example, if ./# m in is slightly lower, as sug gested 
by studies of Type Hp supernova progenitors l lSmartj|2009b . stars 
in this mass range will dominate the total probability due to the 
steepness of the Salpeter IMF. Core masses of the progenitors that 
were not explicitly calculated (./£ < 10 Mq for Z = and 



Jf, < 12 M Q for Z = Zq) can also be different - iNomotol dl984l) 
calculated presupernova structure of a 8.8 Mq star, which has en- 
closed mass at both S/N A - 4fc B and Y c = 0.499 of 1.49 Mq 
(baryonicfl which is signi ficantly higher tha n that assumed for 
low-mass progenitors in the Zhang et al. (2008) models. If we vary 
^min 5 Mq and the baryonic remnant mass of 1.2 M ec ^ 
1.55 Mq for stars where the progenitor structure was not explicitly 
calculated, use uniform priors in ln,/# m i n and lnAf cc , and hold 
Jmax = 100 Mq fixed, we find very little change in the relative 



4 Note that in this model S/N A < 4fc B for all zones below 1.49 Mq 
except a single zone at 1.19 Mq, where it reaches S/N\ = 4.02 fce- 
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probabilities of the models. There i s no significant change of M ec 
with respect to Zha ng et al ] d2008l) . The results on the minimum 
progenitor mass a re typically ^# m in = 8.9 lt{ '\ Mq, whi ch is com- 
patible with both lZhangeTal] d2008h and lSmarttl J2009h . We also 
find 1.32 < Mec < 1. 39 M© for Z = Zq, which is compatible 
with lZhang etail J2008h . 

The current DNS sample has only six systems with precise 
masses. A precise measurement of a NS birth mass greater than 
about 1.5 Mq would greatly constrain the NS formation mod- 
els. Until such system is fountB it might be possible to obtain 
additional constraints by including NS systems with less precise 
mass measurements. The Bayesian formalism naturally accounts 
for non-Ga ussian marginal l ikelihoods Pi(D\ M) such as those pre- 
sented by lOzel et alj (12012) for some systems. Adding DNS sys- 
tems with only a precise total mass measurement is unlikely to 
change our results, as their total masses are compatible with our 
sample (Fig. [5}- We experimented with adding the eclipsing X- 
ray pulsars, which should also have masses near the birth mass 
dRawls et all 1201 it lOzel et al.ll2012h . and found that these mea- 
surements have uncertainties that are too large to improve the con- 
straints. 

The models can also be extended to include other types of 
binaries with degenerate components (e.g. BH-NS, BH-BH, NS- 
WD) since we can calculate the full remnant mass function for the 
binaries (Eq. (8)). That there are no known BH-NS binaries must 
strongly constrain ^# ma x through the relative probabilities of NS- 
NS, NS-BH and BH-BH for different values of .# max . Unfortu- 
nately, this also requires estimates for the relative detection effi- 
ciencies of the individual channels. 

Finally, there are also a number of limitations to this work. The 
current sample of DNSs with precise masses h as only six systems . 
The best available remnant mass function of lZhang etaLM 2008) 
is based on ID models that artificially explode no n-rotating pro- 
genito rs produced by a single stellar evolution code dWooslev et al.l 
120021) . In add ition, specific conditions have to be met to produce a 
DNS system. iBelczvnski et alj d2002l ) gives a comprehensive list 
of possible channels for DNS formation. Most of them involve 
mass transfer and a theoretically uncertain phase of common en- 
velope evolution, which exposes the NS to potentially hypercriti- 



cal accretion (e.g. Ighev alier 1993: iBrownll 1 995t iDewi et al.l F2006 ; 
lLombardi et al.l201ll) . Furthermore. IBelczvnski et al.ld2010h inves- 
tigated the relative numbers of DNS systems and isolated recy- 
cled pulsars and found a disagreement with theoretical predictions 
that point to a lack of understanding of massive binary star evo- 
lution or supernova explosions. Another significant effect is the 
disruption of binaries due to mass ejection and kicks during the 
two supernovae. We implemented the binary survival probability 
after the first sup ernova exp l osion by modifying Equation ((8) us- 
ing the results of Kalogera dl996h . We used the final progenitor 
mass of IZhang et alj d2008l) for the mass of the primary star be- 
fore the explosion and various combinations of the initial and final 
mass for the secondary at the moment of the primary explosion. 
We also tried a number of relative kick velocities. Adding this to 
the calculation had no significant consequence. We note, however, 
that the second supernova is more important for the survival of the 



5 NS masses greater than 1.5 Mq have been meas ured. For examp le, the 
system J1614-2230 has M = 1.97 ± 0.04 M toemorest etal]|2010l) . 
However, these NSs are likely significantly recy cled and do not rep resent 
NS birth masses, although iLin et al.1 J201 ll) and iTauris et al.l J20 1 lh argue 
that the birth mass of J1614-2230 was higher than 1.5 Mq. 



system (e.g. Dewi & van den Heuvell 2004 Willems & Kalogeral 
2004t IWillems et alj l2004t IStairs et al.l l2006t IWang et al.l l2006t 
Wong et al. 2010). Finally, while it is believed the NS masses are 
little affected by mass transfer, this may not hold for the progenitor 
stars. However, we think that results of detailed binary evolution 
and population synthesis models can affect only the details of the 
DNS mass distribution and not the overall conclusions. 



4 DISCUSSIONS & CONCLUSIONS 

In this work, we present a Bayesian framework that directly com- 
pares the observed distribution of NS masses to the theoretical 
models of NS formation in supernovae. We illustrate this method 
on a sample of double neutron stars from lOzel et al.1 d2012h which 
did not experience significant mass accretion and should reflect the 
distribution of NS birth masses. We use the relation between the 
initial progenitor mass and the final remnant mass of Zha ng et"al] 
( 2008), tabulated for a range of explosion energies, primordial and 
solar compositions and for a momentum piston that drives the ex- 
plosion positioned either at the base of the oxygen burning shell 
(the S/Na = 4 &b models) or at the outer edge of the deleptonized 
core (the Y c core models). We also investigate models with no fall- 
back. We assume that the NS progenitors are independently drawn 
from a Salpeter IMF (independent mass model) or from a binary 
model where the primaries (progenitors of millisecond pulsars) are 
drawn from a Salpeter distribution and secondaries (progenitors of 
companions) are drawn from a distribution of mass ratios P(q) that 
is either flat or strongly peaked at q — 1. 

We find a strong preference for binary models over models 
independently drawn from a Salpeter distribution (Fig. |4]. How- 
ever, this is true only if we can correctly assign primary and sec- 
ondary labels to the pulsars and their companions (especially to 
J1906+0746; Fig.[3}. Otherwise, the preference for binary models 
is weak. The mass distributions of primaries and secondaries are 
asymptotically consistent with the IMF, but the mass distribution 
of the primary progenitors is cut off at the minimum mass for NS 
production (./# m i n ), while for secondaries it is cut off at the maxi- 
mum mass (Xax), hence each DNS component probes different 
mean progenitor masses. We do not find any clear preference for 
either of the binary mass ratio distributions we considered. We do 
not find any preference for models with ^# ma x. We find that for 
a primordial composition locating the piston at S/Na = 4 fee is 
strongly favored, while for solar composition models placing the 
piston at the Y c core is more probable. Models with no fallback 
(or with strong explosions that produce little fallback) are always 
preferred because they minimize the width of the NS mass distribu- 
tions (Figs.Q]and[2]l. Globally, the highest probability model has a 
uniform distribution for the binary mass ratios, solar composition, 
sets the remnant mass to that of the Y e core, and has no fallback. 
This model also appears to be consistent with the data in an abso- 
lute "frequentist" sense (Figs. [I] and [5}. 

The preference for no or very little fallback could either be a 
feature of the supernova mechanism or a consequence of mass loss. 
Stripping the progenitor significantly reduces fallbac k, as can be 
seen in the high (^ > 40 Mq) mass solar metallicity Zha ng et "all 
(2008) models. For example, their model with an initial mass of 
60 Mq has pre-explosion mass of only 7.29 M and mass inside 
S/Na = 4 k B of 1.60 Mq, with a fallback mass of only 0.04 M© 
(0.00 M ) for E = 1.2 X 10 51 ergs (2.4 x 10 51 ergs). While the 
mass can be transferred between the stars in the binary and change 
the amount fallback when compared to a single star, the point is 
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Figure 6. Relative probability of the DNS mass models as a function of the 
position of the mass cut in entropy (S/N\ks) that separates the remnant 
and the supernova ejecta (no fallback ). The solid l ine shows the results for 
the supernova progenitors of Wooslev et al. (2002), while the dashed hori- 
zontal line is for the best model based on lZhang et al. ( 2008) that places the 
piston at the Y c core and assumes no fallback (right panel of Fig.|4j- 



that our results prefer no fallback in either star of the binary. This 
is in a greement with the models of IPortegies Zwart & Yungelsonl 
( 1998) where both of the two supernovae in a binary occur in a 
stripped progenitor. Alternatively, the DNS mass distribution may 
represent a surprisingly clear fingerprint of the supernova explo- 
sion mechanism. Specifically, with no fallback, the mass of the NS 
corresponds to the mass coordinate of the progenitor where the 
explosion was initiated. The bulk of the stellar mass growth and 
consequently also the bulk of the DNS population growth occurred 
between redshifts 1 and 2 (e.g. ljuneau et alj20 05), when the metal- 
licity was approximately solar. Our Z — Zq results clearly prefer 
the edge of the deleptonized core (S/Na ~ 2.8 fee) as the mass 
coordinate of the explosion rather than at the base of the oxygen 
shell (S/N A = 4k B ). 

The two piston positions used by IZhang et alj d2008l) arc 
somewhat arbitrary so we explored whether other piston position 
would yie ld better results. We t ook the solar metallicity progenitor 
models of lWooslev et al. (20020 that cover the mass range 10.8 ^ 
^ 40 Mq and determined the mass coordinates that correspond 
to a range of values of the entropy, 2.0 ^ S/Na ^ 4.0 fee. As 
we have no means of calculating the amount of fallback and models 
without fallback are preferred anyway, we assume that the remnant 
masses are set by the mass coordinate at the given entropy. The 
relative probabilities of these models were calculated with ^ m i n 
and M cc as free parameters and ./# ma x fixed (Sec. 13. 3\ and are 
shown as a function of the entropy in Figure [6] We see that there 
is a clear maximum at S/Na ~ 2.8 k B . The relative probability 
at maximum is almost equal to the highest probability model from 
Figure|4] which is not surprising because for many progenitors this 
is the entropy at the edge of the Y e core. If correct, this is an im- 
portant result, because it shows that supernova explosions are initi- 
ated approximately at the edge of the iron core. Since the elapsed 
time from the core bounce until the accretion of all material out 
to the edge of the deleptonized core is 0.1 to 0.3 s for majority of 
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the progenitors, this also means that the supernovae explode by a 
delayed mechanism. That it is not longer challenges some of the 
delayed neutrino-driven mode l s that explode at > 0.5 s after the 
bounce dMarek & Jankdl2009l : llVluller et al.ll2012l) and constrains 
the time a vailable for the un stable modes that potentially drive the 
explosion dFrver et ai1l2012l) . It is important to note that the com- 
position and density of the layers that accrete through the shock 
do not change smoothly with time; the composition interfaces such 
as the one at the edge of the iron core lead to sudden decreases in 
the mass accretion rate, which cause expansion of the shock. This 
expansion mi ght be temporary and can be followed by recession 
of the shock (Marek & Janks i 2009), or it might trigger an explo- 
sion such as those observed after the advection of the Si/SiO shell 
interface dBuras et al.ll2006al;lMuller et al.ll2012l) . 

Independent of any concern about binaries, the greatest limi- 
tation of this and similar studies in the near future is our poor un- 
derstanding of the physics of supernova explosions that link the 
initial progenitor mass to the final remnant mass. While the DNS 
mass distribution strongly suggests that the explosion develops at 
the edge of the deleptonized core, supernova theory does not pro- 
vide any a priori reason why this should be so and whether this is 
true for all progenitor masses, metallicities, rotation rates and other 
parameters. 

To summarize, ideal theoretical calculations of the massive 
star remnant mass function require an understanding of the su- 
pernova explosion mechanism. Until a self-consistent understand- 
ing is available, artificially induced supernova explosions intended 
to study the remnant population should increase their realism by 
manually inducing the explosions by either increasing the neu- 
trino luminosity or absorption cross section in order to prop- 
erly model the transition from accretion to neutrino-driven wind 
and the explosior fl as has been done for some studies of explo- 



m 

tail 



20111; 



2011). Al- 



sion physics (e.g. IScheck et alj|20 06; Nord haus et 
iFuiimoto et alj|201 ll ; Irlanke et alj|201 ll ; iKotake et 
though the DNS distribution suggests that fallback cannot be sig- 
nificant for the DNS systems, it would be interesting to see if this 
is also a feature of "semi-physical" explosion models. Of particu- 
lar importance is whether fallback is a strong or stochastic func- 
tion of progenitor structure, if it occurs. The range of progenitor 
models and model supernovae also needs to be expanded to better 
cover the progenitor mass range in order to fully sample the rapidly 
changing core properties with mass. Models are most needed near 
solar metallicities, since such stars overwhelmingly dominate any 
observable population of supernovae or NS binaries, and at the low 
masses that dominate the progenitor population due to the initial 
mass function. Observations mainly constrain the outcomes of low- 
mass solar metallicity stars - the models least examined in theoret- 
ical studies. 
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7 After the submission of this paper, Ugliano et al. ( 2012) published a NS 
mass distribution that is based on neutrino-driven explosions that are fol- 
lowed until the end of the fallback. We present the resulting probability 
distribution and comparison with other models in the Appendix. 
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Table 1. Summary of the remnant mass distribution models from Zhang et al. (2008). 





l iston cit 


•*min 


Progenitor mass range 


rp [1 r\51 qvo-q! 
XV [1U <JI£>oJ 





S/N A =4k B 


9.5 Mq 


(10, 100) Mq 


0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.4, 3.0, 5.0, 10.0 





S/N A = 4k B 


9.5 M 


(10, 100) Mq 


S/N A = 4 fcg core mass, no fallback 





Y e core 


9.5 M 


(10, 100) M 


1.2, 10.0 





Y c core 


9.5 Mq 


(10, 100) Mq 


Y e core mass, no fallback 


z e 


S/N A = 4fe B 


9.1 Mq 


(12, 100) Mq 


1.2, 2.4 


z & 


S/N A = 4fc B 


9.1 Mq 


(12, 100) Mq 


S/N A = 4 &b core mass, no fallback 


% 


Y core 


9.1 Mq 


(12, 100) Mq 


1.2, 2.4 


z & 


Y c core 


9.1 Mq 


(12, 100) Mq 


Y e core mass, no fallback 



Table 2. Summary of double neutron star systems. 



Name 


Mi 




Pi [ms] 


Pi 




M 2 


C2 


J0737-3039 a 


1.3381 


0.0007 


22.7 


1.8 x 10" 


18 


1.2489 


0.0007 


B1534+12 


1.3332 


0.0010 


37.9 


2.4 x 10" 


18 


1.3452 


0.0010 


J1756-2251 


1.40 


0.02 


28.5 


1.0 x 10" 


18 


1.18 


0.02 


J 1906+0746 


1.248 


0.018 


144.1 


2.0 x 10" 


14 


1.365 


0.018 


B1913+16 


1.4398 


0.002 


59.0 


8.6 x 10" 


18 


1.3886 


0.002 


B2127+11C 


1.358 


0.010 


30.5 


5.0 x 10" 


18 


1.354 


0.010 
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a The companion is also a pulsar with P2 = 2.8 s and P2 
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APPENDIX A: THE NS MASS DISTRIBUTION OF 
UGLIANO ET AL. (2012) 

After the submission of this paper, Uglian o et al. I d2012h published 
a NS mass distribution that is based on ID simulations of neutrino- 
driven explosions followed from the onset of the collapse until the 
end of the fallback. The simulations are normalized by compari- 
son with the observed parameters of SN 1987A. The distribution is 
based on over 100 progenitors with 10 M© ^ ^£ ^ 40 Mq, which 
we supplement with remnant masses of 1.35 Mq for 9.1 M q ^ 
J( < 10 Mq, similarly to the models of Izhang etH\ d2008h dis- 
cussed in the main paper. We also corrected the remnant masses 
for the binding energy according to Equation l(9)- In Figure [ATI we 
show the resulting probability distribution of the DNS masses. We 
see that the peaks of the distribution are shifted by ~ 0.1 Mq to 
higher masses with the respect to the DNS data. Even though the 
fallback on the remnants is included, the probability distribution 
has little power for NS masses higher than ab out 1.6 Mq, unlike 
the large amounts of fallb ack in many of the Izhang et all (2008) 
models (Fig. 0. Indeed, lUgliano et all d2012h find little fallback 
for high-mass progenitors. 

In order to compare the relative probability of the 
lUgliano et al.1 d2012h distribution to the models in Figure [6] we 
marginalize over ^ m i n and the remnant mass M cc for progeni- 
tors with ^# min ^ ^ < 10.8 Mq in the same way as we did for 
the models of lWooslev etai] ( l2002h . We find that the relative prob- 
ability is only 10 3 ' 2 on the scale of Figure|6] which is likely caused 
by the slight mismatch in the probability peaks with respect to the 
DNS data. If we include remnant masses without fallback, which 
might be more appropriate for DNS progenitors with stripped hy- 
drogen envelopes, the relative probability raises to about 10 5 ' 4 . 
However, this is about a factor of 10 worse than our best models in 
Figure s [4]and[6] Thus, the reduction of fallback in thelU gliano et al.l 
( 20 1 2) NS mass distribution provides a better match to the observed 
masse s of DNS binaries, but not as good match as the Zha ng et al.l 
d2008h Y c core models with no fallback. 



© 0000 RAS, MNRAS 000, 000-000 



14 Pejcha, Thompson and Kochanek 




millisecond pulsar mass M A [M Q ] ?(¥jll ma J 

Figure Al. Same as Figure[T] but for a model based on the solar m etallicity NS mass distribution of lUgliano et aljfeplA . The distribution of P(g) is uniform 
and the calculation includes all progenitors o f lUglianoetaflfeOl^ . 
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